Wczytanie zbioru danych Mammographic Mass

Zbiór danych wczytany został do tablicy data. Niekompletne dane zostały zamienione na wartość NaN, a pozostałe - na wartości numeric.

#wczytanie danych
data<-read.table("mammographic_masses.data", header=FALSE, sep=",");
colnames(data) <- c('bi_rads', 'age', 'shape', 'margin', 'density', 'severity')


nan_data <- which(data=='?');
new_data <- replace(data, data=='?', NA); #zamiana ? na NA
new_data <- as.data.frame(lapply(new_data,as.numeric)) # konwersja wszystkich wartosci w dataframie na numeric

Analiza zbioru danych Mammographic Mass

Histogramy dla poszczególnych cech w zbiorze danych

Na podstawie poniższych histogramów można dokonać analizy statystycznej cech, np. wyznaczyć wartości minimalne i maksymalne oraz zauważyć wartości odstające.
W celu ujednolicenia wykresów, do konkretnych cech przypisane zostały kolory:

# Bindowanie kolorow do cech
featureColors <- c(bi_rads="#9ACD32",
                  age="#F08080",
                  shape="#B0E0E6",
                  margin="#DEB887",
                  density="#BA55D3",
                  severity="#66CDAA")


Histogram cechy BI-RADS

bi_rads_h <- hist(new_data$bi_rads, main="Ocena BI-RADS", xlab="Wartosc", ylab="Ilosc wystapien", xlim=c(0,max(new_data$bi_rads, na.rm = TRUE)), ylim=c(0, 600), breaks=56, col=featureColors['bi_rads'])

zero = which(bi_rads_h$counts == 0)
text(bi_rads_h$mids[-zero],bi_rads_h$counts[-zero],labels=bi_rads_h$counts[-zero], adj=c(0.5, -0.5))


Histogram cechy PATIENT’S AGE

age_h <- hist(new_data$age, main="Wiek pacjentki (w latach)", xlab="Wiek", ylab="Ilosc wystapien", xlim=c(0,100), breaks=99, col=featureColors['age'])


Histogram cechy MASS SHAPE

shape_h <- hist(new_data$shape, main="Ksztalt guza", xlab="Ksztalt", ylab="Ilosc wystapien", xlim=c(0,4), ylim=c(0, 450), breaks=0:4, col=featureColors['shape'], xaxt="n")
text(shape_h$mids,shape_h$counts,labels=shape_h$counts, adj=c(0.5, -0.5))
axis(1, shape_h$mids, c("round", "oval", "lobular", "irregular"), line=-1, lwd=0, lwd.ticks = 1)


Histogram cechy MASS MARGIN

margin_h <- hist(new_data$margin, main="Margines guza", xlab="Margines", ylab="Ilosc wystapien", xlim=c(0,max(new_data$margin, na.rm = TRUE)), ylim=c(0, 400), breaks=0:5, col=featureColors['margin'], xaxt="n")
text(margin_h$mids,margin_h$counts,labels=margin_h$counts, adj=c(0.5, -0.5))
axis(1, at=margin_h$mids, labels = rep_len("", length(margin_h$mids)), lwd=0, lwd.ticks = 1, line = -1)
text(margin_h$mids, par("usr")[3] - 0.2, labels = c("circumscribed", "microlobulated", "obscured", "ill-defined", "spiculated"),
     srt = 35, pos = 1, xpd = TRUE)


Histogram cechy MASS DENSITY

density_h <- hist(new_data$density, main="Gestosc guza", xlab="Gestosc guza", ylab="Ilosc wystapien", xlim=c(0,max(new_data$density, na.rm = TRUE)), ylim=c(0, 850), breaks=0:4, col=featureColors['density'], xaxt="n")
text(density_h$mids,density_h$counts,labels=density_h$counts, adj=c(0.5, -0.5))
axis(1, density_h$mids, c("high", "iso", "low", "fat-containing"), line=-1, lwd=0, lwd.ticks = 1)


Histogram cechy MASS SEVERITY

severity_h <- hist(new_data$severity, right=F, main="Stopien zaawansowania guza", xlab="Zlosliwosc", ylab="Ilosc wystapien", xlim=c(0,1), ylim = c(0,600), breaks=2, col=featureColors['severity'], xaxt="n")
text(severity_h$mids,severity_h$counts,labels=severity_h$counts, adj=c(0.5, -0.5))
axis(1, severity_h$mids, c("benign", "malignant"), line=-1, lwd=0, lwd.ticks = 1)

Statystyki zbioru danych Mammographic Mass (z pominieciem wartości NA)

Poniżej zaimplementowana funkcja static_values służy do wyznaczania wartości statystycznych poszczególnych cech, z pominięciem wartości NaN. Oblicza ona wartości minimalne i maksymalne, średnią, odchylenie standardowe i wariancję, medianę oraz dominantę.

getmode <- function(v, na.rm=TRUE) {
  uniqv <- unique(v)
  uniqv[which.max(tabulate(match(v, uniqv)))]
}

static_values <- function(v) {
  min <- min(v, na.rm = TRUE)
  max <- max(v, na.rm = TRUE)
  mean <- mean(v, na.rm = TRUE)
  sd <- sd(v, na.rm = TRUE)
  var <- var(v, na.rm = TRUE)
  median <- median(v, na.rm = TRUE)
  mode <- getmode(v, na.rm = TRUE)
  res <- c(min, max, mean, median, mode, sd, var)
  names(res) <-c('Min', 'Max', 'Mean', 'Median', 'Mode', 'Sd', 'Var')
  return(res)
}
#--------------------------------------------------
#wartosci statystyczne zbioru danych
(data_stats <- sapply(new_data, static_values))
##          bi_rads       age    shape   margin   density  severity
## Min     0.000000  18.00000 1.000000 1.000000 1.0000000 0.0000000
## Max    55.000000  96.00000 4.000000 5.000000 4.0000000 1.0000000
## Mean    4.348279  55.48745 2.721505 2.796276 2.9107345 0.4630593
## Median  4.000000  57.00000 3.000000 3.000000 3.0000000 0.0000000
## Mode    4.000000  59.00000 4.000000 1.000000 3.0000000 0.0000000
## Sd      1.783031  14.48013 1.242792 1.566546 0.3804439 0.4988932
## Var     3.179201 209.67419 1.544532 2.454065 0.1447376 0.2488944

Porządkowanie zbioru danych

Porządkowanie zbioru danych miało na celu usunięcię niekompletnych danych oraz określenie postępowania w przypadku wartości odstających.

Zbiór danych z usuniętymi wartościami NaN zapisany został w zmiennej data_clean.

Jedyną zauważoną wartością odstającą jest wartość 55 określająca cechę BI-RADS. Prawdopodobnie powstała ona w wyniku błędu przy przepisywaniu danych. Można założyć, że liczba 55 jest zdupliowaną cyfrą 5. Z tego powodu w zbiorze Mammographic Mass została one ostatecznie zapisana jako wartość 5. Wpływ na to miał równiez fakt, że wektor z ta wartością nie zawierał żadnych niekompletnych danych, a po wyczyszczeniu zbioru jego liczebność zmalała do 830. usunięcie kompletnego wektoru mogłoby narazić sieć na utratę istotnych danych.

#usunięcie niekompletnych danych
data_clean <- new_data[complete.cases(new_data), ]


#element odstajacy
data_clean$bi_rads[data_clean$bi_rads==55] <- 5

Analiza zbioru danych Mammographic Mass po wyczyszczeniu danych

Histogramy dla poszczególnych cech w wyczyszczonym zbiorze danych


Histogram cechy BI-RADS

bi_rads_new_h <- hist(data_clean$bi_rads, main="Ocena BI-RADS", xlab="Wartosc", ylab="Ilosc wystapien", xlim=c(0,6), ylim=c(0, 500), breaks=-0.5:6.5, col=featureColors['bi_rads'], xaxt="n")
text(bi_rads_new_h$mids,bi_rads_new_h$counts,labels=bi_rads_new_h$counts, adj=c(0.5, -0.5))
axis(1, bi_rads_new_h$mids, seq_along(bi_rads_new_h$mids))


Histogram cechy PATIENT’S AGE

age_new_h <- hist(data_clean$age, main="Wiek pacjentki (w latach)", xlab="Wiek", ylab="Ilosc wystapien", xlim=c(0,100), breaks=99, col=featureColors['age'])


Histogram cechy MASS SHAPE

shape_new_h <- hist(data_clean$shape, main="Ksztalt guza", xlab="Ksztalt", ylab="Ilosc wystapien", xlim=c(0,4), ylim=c(0,400), breaks=0:4, col=featureColors['shape'], xaxt="n")
text(shape_new_h$mids,shape_new_h$counts,labels=shape_new_h$counts, adj=c(0.5, -0.5))
axis(1, shape_new_h$mids, c("round", "oval", "lobular", "irregular"), line=-1, lwd=0, lwd.ticks = 1)


Histogram cechy MASS MARGIN

margin_new_h <- hist(data_clean$margin, main="Margines guza", xlab="Margines", ylab="Ilosc wystapien", xlim=c(0,max(data_clean$margin, na.rm = TRUE)), ylim=c(0,400), breaks=0:5, col=featureColors['margin'], xaxt="n")
text(margin_new_h$mids,margin_new_h$counts,labels=margin_new_h$counts, adj=c(0.5, -0.5))
axis(1, at=margin_new_h$mids, labels = rep_len("", length(margin_new_h$mids)), lwd=0, lwd.ticks = 1, line = -1)
text(margin_new_h$mids, par("usr")[3] - 0.2, labels = c("circumscribed", "microlobulated", "obscured", "ill-defined", "spiculated"),
     srt = 25, pos = 1, xpd = TRUE)


Histogram cechy MASS DENSITY

density_new_h <- hist(data_clean$density, main="Gestosc guza", xlab="Gestosc guza", ylab="Ilosc wystapien", xlim=c(0,max(data_clean$density, na.rm = TRUE)), ylim=c(0,800), breaks=0:4, col=featureColors['density'], xaxt="n")
text(density_new_h$mids,density_new_h$counts,labels=density_new_h$counts, adj=c(0.5, -0.5))
axis(1, density_h$mids, c("high", "iso", "low", "fat-containing"), line=-1, lwd=0, lwd.ticks = 1)


Histogram cechy MASS SEVERITY

severity_new_h <- hist(data_clean$severity, right=F, main="Stopien zaawansowania guza", xlab="Zlosliwosc", ylab="Ilosc wystapien", xlim=c(0,1), ylim = c(0,500), breaks=2, col=featureColors['severity'], xaxt="n")
text(severity_new_h$mids,severity_new_h$counts,labels=severity_new_h$counts, adj=c(0.5, -0.5))
axis(1, severity_new_h$mids, c("benign", "malignant"), line=-1, lwd=0, lwd.ticks = 1)

Analiza zbioru pod katem determinant zlosliwosci guza

Sprawdzono udział guzów łagodnych i złośliwych w zbiorze, w zależności od wartości danej cechy i zwizualizowano w postaci barplotów. Wgląd w wykreślone dane pozwala (w pewnym ograniczonym stopniu) przewidywać, dla których wartości których cech, spodziewać się możemy, że masa guzowa będzie łagodna, a która złośliwa.

bi_rads_benVSmag <- as.matrix(t(prop.table(table(data_clean[c('bi_rads', 'severity')]), 1)))
bi_rads_benVSmag <- cbind(bi_rads_benVSmag[,1], c(0.0,0.0), bi_rads_benVSmag[,-1])
colnames(bi_rads_benVSmag) <- as.character(0:6)

age_benVSmag <- t(prop.table(table(data_clean[c('age', 'severity')]), 1))
shape_benVSmag <- t(prop.table(table(data_clean[c('shape', 'severity')]), 1))
margin_benVSmag <- t(prop.table(table(data_clean[c('margin', 'severity')]), 1))
density_benVSmag <- t(prop.table(table(data_clean[c('density', 'severity')]), 1))


# barploty
par(mar = c(5.1, 4.1, 4.1, 6))
b <- barplot(bi_rads_benVSmag, col = c("seagreen3", "firebrick3"),
             main = "Ocena BI-RADS vs. udzial guzow lagodnych/zlosliwych",
             xlab = "Ocena BI-RADS", ylab = "%",
             legend.text = c("benign", "malignant"),
             args.legend = list(x = "topright", inset = c(-0.2, 0)))
text(b[2], 0.48, "BRAK", col="tomato2", srt=90, font=2)

barplot(age_benVSmag, col = c("seagreen3", "firebrick3"),
        main = "Wiek pacjentki vs. udzial guzów lagodnych/zlosliwych",
        xlab = "Wiek pacjentki", ylab = "%",
        legend.text = c("benign", "malignant"),
        args.legend = list(x = "topright", inset = c(-0.2, 0)))

b <- barplot(shape_benVSmag, col = c("seagreen3", "firebrick3"),
             main = "Ksztalt masy vs. udzial guzow lagodnych/zlosliwych",
             xlab = "Ksztalt masy", ylab = "%", xaxt="n",
             legend.text = c("benign", "malignant"),
             args.legend = list(x = "topright", inset = c(-0.2, 0)))

axis(1, b, c("round", "oval", "lobular", "irregular"), line=-1, lwd=0, lwd.ticks = 1)

b <- barplot(margin_benVSmag, col = c("seagreen3", "firebrick3"),
             main = "Margines masy vs. udzial guzow lagodnych/zlosliwych",
             xlab = "Margines masy", ylab = "%", xaxt="n",
             legend.text = c("benign", "malignant"),
             args.legend = list(x = "topright", inset = c(-0.2, 0)))

axis(1, b, rep("", length(b)), line=-1, lwd=0, lwd.ticks = 1)
text(b, par("usr")[3], labels = c("circumscribed", "microlobulated", "obscured", "ill-defined", "spiculated"),
     srt = 25, pos = 1, xpd = TRUE)

b <- barplot(density_benVSmag, col = c("seagreen3", "firebrick3"),
             main = "Gestosc masy vs. udzial guzow lagodnych/zlosliwych",
             xlab = "Gestosc masy", ylab = "%", xaxt="n",
             legend.text = c("benign", "malignant"),
             args.legend = list(x = "topright", inset = c(-0.2, 0)))

axis(1, b, c("high", "iso", "low", "fat-containing"), line=-1, lwd=0, lwd.ticks = 1)

Macierz korelacji

Macierz korelacji pozwala zwizualizować korelacje między cechami zbioru danych. Jako miarę korelacji między cechami zastosowano współczynnik korelacji Pearsona. Z wykresu usunięto elementy znajdujące się ponad główną przekątną macierzy, aby niepotrzebnie nie duplikować korelacji między poszczególnymi cechami.

require(reshape2)
require(plotly)
cormat <- round(cor(data_clean), 2);
melted_cormat <- melt(cormat);

# Get lower triangle of the correlation matrix
get_lower_tri<-function(cormat){
  cormat[upper.tri(cormat)] <- NA
  return(cormat)
}
# Get upper triangle of the correlation matrix
get_upper_tri <- function(cormat){
  cormat[lower.tri(cormat)]<- NA
  return(cormat)
}

upper_tri <- get_upper_tri(cormat);

reorder_cormat <- function(cormat){
  # Use correlation between variables as distance
  dd <- as.dist((1-cormat)/2)
  hc <- hclust(dd)
  cormat <-cormat[hc$order, hc$order]
}

cormat <- reorder_cormat(cormat)
upper_tri <- get_upper_tri(cormat)
# Melt the correlation matrix
melted_cormat <- melt(upper_tri, na.rm = TRUE)

ggheatmap <- ggplot(melted_cormat, aes(Var2, Var1, fill = value))+
      geom_tile(color = "white")+
      scale_fill_gradient2(low = "black", high = "firebrick2", mid = "dodgerblue4", 
                           midpoint = 0, limit = c(-1,1), space = "Lab", 
                           name="Pearson\nCorrelation") +
      theme_minimal()+ # minimal theme
      theme(axis.text.x = element_text(angle = 45, vjust = 1, 
                                       size = 12, hjust = 1))+
      theme(axis.text.y = element_text(size = 12)) + 
      coord_fixed()
    
    ggplotly(ggheatmap + 
      geom_text(aes(Var2, Var1, label = value), color = "white", size = 4) +
      theme(
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank(),
        axis.ticks = element_blank(),
        legend.justification = c(1, 0),
        legend.position = c(0.6, 0.7),
        legend.direction = "horizontal")+
      guides(fill = guide_colorbar(barwidth = 7, barheight = 1,
                                   title.position = "top", title.hjust = 0.5)))

Standaryzacja zbioru danych

Do standardowego rozkładu normalnego N(0, 1), bez labeli (czyli bez cechy severity). Standardowy rozkład normalny każdej cechy uzyskujemy poprzez odjęcie od wartości każdej cechy w zbiorze jej średnią wartość, a następnie podzielenie uzyskanego wyniku przez odchylenie standardowe cechy. Wyznaczone średnie (scaled_centers) i odchylenia standardowe (scaled_scales) cech zapamiętujemy w zmiennych w celu standaryzowania za ich pomocą danych dostarczonych przez użytkownika.

s <- scale(data_clean[,2:5])

scaled_centers <- attr(s, 'scaled:center')
scaled_scales <- attr(s, 'scaled:scale')

data_stand <- data.frame(age=s[,1], shape=s[,2], margin=s[,3], density=s[,4], data_clean['severity'])
head(data_stand, 10)
##            age      shape     margin    density severity
## 1   0.76460188  0.1755305  1.3953435  0.2403212        1
## 3   0.15117947  0.9804496  1.3953435  0.2403212        1
## 4  -1.89356190 -1.4343076 -1.1570204  0.2403212        0
## 9   0.08302143 -1.4343076  1.3953435  0.2403212        1
## 11  1.37802429 -1.4343076  0.7572525  0.2403212        1
## 12 -0.93934926 -0.6293885 -1.1570204  0.2403212        1
## 14 -1.34829753  0.1755305 -1.1570204 -2.6092013        0
## 15  0.28749556 -0.6293885 -1.1570204 -2.6092013        0
## 16 -0.12145271 -1.4343076 -1.1570204  0.2403212        0
## 17 -0.25776880  0.1755305  0.7572525  0.2403212        0

Utworzenie oraz trening sieci neuronowej

Tworzymy sieć złożoną z 1 warstwy ukrytej o 6 neuronach. Zastosowaną funkcją aktywacji jest funkcja sigmoid unipolarna, natomiast zastosowaną funkcją błędu jest błąd średniokwadratowy. Zbiór ustandaryzowanych danych podzielony został na zbiory treningowy oraz testowy w stosunku 0.8:0.2. Trening sieci dokonywany był na zbiorze treningowym. Podział danych z wykorzystaniem funkcji createDataPartition z pakietu caret wykonywany jest z zachowaniem proporcji klas z zbiorach (podział ze stratyfikacją). Jako klasę “pozytywną” przyjęto klasę ‘malignant’.

library(rpart)
library(caret)
require(neuralnet)

#podział na zbior treningnowy i testowy
#sev <- unlist(data_clean['severity'])


sev <- data_clean$severity
set.seed(100)
part <- createDataPartition(sev, p = 0.8, list  = FALSE)

training_set <- scale(data_clean[part, 2:5])
scaled_centers <- attr(training_set, 'scaled:center')
scaled_scales <- attr(training_set, 'scaled:scale')

testing_set <- data_clean[-part, 2:5]

training_sev <-sev[part]
testing_sev <-sev[-part]



data_stand_tr <- data.frame(age=training_set[,1], shape=training_set[,2], margin=training_set[,3], density=training_set[,4], training_sev)
data_ts <- data.frame(age=testing_set[,1], shape=testing_set[,2], margin=testing_set[,3], density=testing_set[,4], testing_sev)

nn <- neuralnet(training_sev ~ age+shape+margin+density, data=data_stand_tr,
                err.fct = "sse", hidden = c(6), act.fct = "logistic")

predsVStarget_tr <- data.frame(case=rownames(data_stand_tr),
                               predictions=factor(round(unlist(nn$net.result), digits = 0), labels = c('benign', 'malignant')),
                               target=factor(data_stand_tr$training_sev, labels = c('benign', 'malignant')))

pred <-round(predict(nn, scale(testing_set, center = scaled_centers, scale = scaled_scales)), 2)

predsVStarget_ts <-data.frame(case = rownames(testing_set), 
                              predictions=factor(round(pred, digits = 0), labels = c('benign', 'malignant')),
                              target = factor(data_ts$testing_sev, labels = c('benign', 'malignant')))

confMatrix_tr <- confusionMatrix(data=predsVStarget_tr$predictions, reference=predsVStarget_tr$target, positive = "malignant")
confMatrix_ts <- confusionMatrix(data=predsVStarget_ts$predictions, reference=predsVStarget_ts$target, positive = "malignant")

Macierze pomyłek

Macierze pomyłek służą do oceny osiągów sieci zarówno na zbiorze danych treningowych jak i zbiorze danych testowych. Zawierają one m.in. informacje o dokładności predykcji sieci, a także ocenę osiągniętej czułości oraz specyficzności.

confMatrix_tr
## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  benign malignant
##   benign       276        47
##   malignant     69       272
##                                           
##                Accuracy : 0.8253          
##                  95% CI : (0.7942, 0.8534)
##     No Information Rate : 0.5196          
##     P-Value [Acc > NIR] : <2e-16          
##                                           
##                   Kappa : 0.651           
##                                           
##  Mcnemar's Test P-Value : 0.0512          
##                                           
##             Sensitivity : 0.8527          
##             Specificity : 0.8000          
##          Pos Pred Value : 0.7977          
##          Neg Pred Value : 0.8545          
##              Prevalence : 0.4804          
##          Detection Rate : 0.4096          
##    Detection Prevalence : 0.5136          
##       Balanced Accuracy : 0.8263          
##                                           
##        'Positive' Class : malignant       
## 
fourfoldplot(confMatrix_tr$table, color = c("firebrick2", "seagreen3"),
             conf.level = 0, margin = 1, main = "Macierz pomylek na zbiorze trenigowym")

confMatrix_ts
## Confusion Matrix and Statistics
## 
##            Reference
## Prediction  benign malignant
##   benign        60         9
##   malignant     22        75
##                                           
##                Accuracy : 0.8133          
##                  95% CI : (0.7455, 0.8694)
##     No Information Rate : 0.506           
##     P-Value [Acc > NIR] : < 2e-16         
##                                           
##                   Kappa : 0.6257          
##                                           
##  Mcnemar's Test P-Value : 0.03114         
##                                           
##             Sensitivity : 0.8929          
##             Specificity : 0.7317          
##          Pos Pred Value : 0.7732          
##          Neg Pred Value : 0.8696          
##              Prevalence : 0.5060          
##          Detection Rate : 0.4518          
##    Detection Prevalence : 0.5843          
##       Balanced Accuracy : 0.8123          
##                                           
##        'Positive' Class : malignant       
## 
fourfoldplot(confMatrix_ts$table, color = c("firebrick2", "seagreen3"),
             conf.level = 0, margin = 1, main = "Macierz pomylek na zbiorze testowym")